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ABSTRACT 

One of the most interesting and challenging aspects of formation guidance law design is the coupling of the orbit 
design and the science return. The analyst’s role is more complicated than simply to design the formation geometry 
and evolution. He or she is also involved in designing a significant portion of the science instrument itself. The 
effectiveness of the formation as a science instrument is intimately coupled with the relative geometry and evolution of 
the collection of spacecraft. Therefore, the science return can be maximized by optimizing the orbit design according 
to a performance metric relevant to the science mission goals. In this work, we present a simple method for optimal 
formation guidance that is applicable to missions whose performance metric, requirements, and constraints can be 
cast as functions that are explicitly dependent upon the orbit states and spacecraft relative positions and velocities. 
We present a general form for the cost and constraint functions, and derive their semi-analytic gradients with respect 
to the formation initial conditions. The gradients are broken down into two types. The first type are gradients of 
the mission specific performance metric with respect to formation geometry. The second type are derivatives of the 
formation geometry with respect to the orbit initial conditions. The fact that these two types of derivatives appear 
separately allows us to derive and implement a general framework that requires minimal modification to be applied to 
different missions or mission phases. To illustrate the applicability of the approach, we conclude with applications to 
two missions: the Magnetospheric Multiscale mission (MMS), and the Laser Interferometer Space Antenna (LISA). 

INTRODUCTION 

Developing guidance law algorithms for spacecraft formations and constellations is a challenging problem that can 
seemingly involve mutually exclusive goals. From the perspective of a mission analyst working on orbit design for 
distributed spacecraft missions, it is desirable to have a flexible method that can be applied to many problems, and 
that permits modification of the cost and constraints with as little effort as possible. We also desire an approach that 
can handle real-world mission constraints, and an optimal solution in a reasonable amount of time, and converge 
consistently with relatively poor initial guesses. The work presented in this paper was motivated by the need for 
an algorithm that can provide formation initial conditions for many different phases in the mission design process. 
These phases include, but are not limited to: developing initial conditions to use as target states to achieve after 
separation from the launch vehicle and separation from the spacecraft stack, developing new conditions for formation 
maintenance after perturbations and errors in navigation and control have caused a significant degradation in mission 
return, and developing new formation initial conditions for a new mission phase when the desired formation geometry 
is significantly different. For each of these mission phases we need to determine initial conditions for the orbits of 
multiple, cooperating spacecraft, that optimally satisfy mission requirements and constraints. This is a long list of 
goals, but we can come very close to meeting them all by developing effective design algorithms. 

The approach presented in this paper is a direct parameter optimization method which assumes that at a given 
instant in time, one can formulate a performance metric that is an explicit function of the inertial cartesian states, 
and the relative position vectors, relative velocity vectors, range and range rates of all spacecraft in formation. The 
performance metric allows us to quantify the effectiveness of the relative dynamics of the spacecraft at a particular 
point in the orbit. The cost function, J, is simply the integral of the instantaneous performance metric over regions 
of interest to the particular mission. The goal is to maximize or minimize the cost function, according to the specific 
mission, while simultaneously satisfying equality and inequality constraints consistent with mission requirements and 
real-world limitations. 

We begin by posing a continuous-time set of cost and constraint functions that define a general framework for the 
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problem. The continuous-time system is the problem we would like to solve. However, there are certain advantages 
to discretizing the problem. Hence, we discretize the problem to enable a Nonlinear Programming (NLP) routine 
to solve for an optimal solution. The discrete cost function can be calculated by quadrature, or a simple trapezoid 
rule, using numerical integration of the non-linear orbit equations of motion, with perturbations, to give us the 
state of the formation at discrete times in the orbit. The constraints are calculated in a similar manner. We show 
how to calculate the derivative of the cost function and the Jacobian of the constraint functions using the State 
Transition Matrix (STM). The STM is found by numerically integrating an additional 36 equations per spacecraft 
with terms that include perturbations. There is a significant advantage to this approach. Numerical integration, 
while non-trivial, is a solved problem and something that modern computers perform efficiently. Secondly, as we 
will see in a later section, the derivatives of the cost and constraints contain two types of terms. The first type of 
terms are composed of portions of the STM and inertial states. The second type of terms contain derivatives of 
algebraic functions with respect to the formation geometry. Because the terms are separated in this way, we can 
implement the algorithm in a very general way so as to allow minimal modification to new problems. We conclude 
the paper with two applications that illustrate how the method is applied to real-world problems. The first example 
is the Magnetospheric Multi- Scale (MMS) mission. The second example is the Laser Interferometer Space Antenna 
(LISA). Let’s begin by defining some notation. 

SYMBOLS 

Variables 

r Position vector 

v Velocity vector 

n Number of spacecraft 

m Number of unique sides in the formation 

N Number of path constraints 

rik Number of quadrature points 

Cfc Quadrature constant at point k 

s ik Vector defining side i at quadrature point k 

s ik Rate vector of side z, at point k 

Sik Length of side z, at point k 

Sik Rate of change side z, at point k 

$ Orbit state transition matrix 

A Upper left 3x3 partition of & 

B Upper right 3x3 partition of $ 

C Lower left 3x3 partition of $ 

D Lower right 3x3 partition of $ 

Subscripts 

i Side index 

k Quadrature point index 

j Spacecraft index 

£ Spacecraft index 

p Dummy index 

J Indicates association with cost function 

C Indicates association with constraints 

General Problem Statement 

The reliability and effectiveness of an optimization algorithm is intimately dependent upon the parameterization of 
the cost and constraint functions. In this section we pose a set of cost and constraint functions for a formation of 
spacecraft as explicit functions of the cartesian states of the spacecraft in formation, and the relative motion of the 
spacecraft. Let’s begin with some definitions. Recall that our goal is to find a set of initial cartesian states for a set 
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of n spacecraft. The vector defining the initial cartesian states for all spacecraft is 


= 

■ol 
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where 


Xoi = [rji v^i] T = [x 0 i Vol Zol ®ol jfol *olP 


( 2 ) 


and so on. For our implementation, we chose to use the Earth MJ2000 equatorial coordinate system to express the 
cartesian states. However, it is possible to work in other coordinate systems, including rotating systems, as long as 
you are consistent. We assume that the equations of motion for the spacecraft are explicit functions of the spacecraft 
state and time, or: 

X(t) = f(X,t) (3) 

Let’s define the i th unique side vector of the formation, s*, as the vector from the £ th to the j th spacecraft 


s i(t) = r j{t) -r e {t) 


(4) 


The subscripts j and i can be chosen arbitrarily as long as for each value of i there is a given j y £ pair, the pairs 
result in a unique side, 1 < j < n and 1 < £ < n, and the chosen definition is used consistently throughout the 
implementation. Similarly to the side vector, the side rate vector s n the side length s*, and the side rate Si are 
defined as 


s,(t) =r j(t) -r e(t) 
Si{t) = (si(t) T Si(t )) 1/2 


(5) 

(6) 
(7) 


Let’s pose a continuous-time cost function of the form 



Fj (X(t), Si ( t ) , Si(t) , Si ( t ), Si (t) , C j) dt 


( 8 ) 


where C j is a set of constants. Note that we assume the cost function is an explicit function of the relative geometry 
of the spacecraft. In Eq. (8), Fj is the algebraic performance metric that contains information on how well the 
geometry satisfies mission goals at a particular instant. Fj is mission specific and must be formulated according 
to a given problem. In a later section we present two examples for illustrative purposes. Similarly, we can define 
nonlinear equality constraints G#, and nonlinear inequality constraints, G/, as 


G e(X 0 ) =FE(X(t),s i (t),s i (t) i s i (t),Si(t),c E ) = C e (9) 

G/(X 0 ) = Fj (X(t), Si(t),Si(t)j Si(t), Si(t), cj) <Cj (10) 

Finally, we also impose linear equality constraints, linear inequality constraints, and bound constraints, respectively 
as 

A eX c = (11) 

A/X 0 < Bj (12) 

L < X 0 < U (13) 


We seek a solution that minimizes J, while simultaneously satisfying the constraints in Eqs. (9-13). However, 
there are several issues that make this difficult, if not impossible, to solve exactly. The first is that it is not possible 
to develop closed form solutions of X in terms of t. Similary, we don’t know how to express the relative geometry 
as explicit functions of time. Certain nonlinear constraint functions can also cause difficulties. By discretizing 
the problem, we can avoid these complications. Yet, we must be satisfied with a solution that is not exactly 
optimal. However, the difference between the rigorously optimal solution and the near-optimal solution obtained 
from parameter optimization is often so small as to make the difference academic. Now let’s take a look at how to 
discretize the continuous-time problem. 
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Discretization of the Cost and Constraints 


In general, we cannot solve the integral in Eq. (8). Also, it is often difficult to evaluate the path constraint function 
in Eqs. (9-13) along the entire trajectory. Furthermore, differentiating the equations with respect to X G is usually 
difficult. For these reasons, it is convenient to approximate the integral in the cost function by using a quadrature 
rule. Similarly, it is convenient to evaluate path constraints at discrete points along the trajectory. Once we have 
recast the cost and constraints as a discrete system, we can take the derivatives relatively easy as we’ll see near the 
end of this section. 


Let’s begin by writing the cost function as a quadrature rule as follows: 

n-k 

»/(X 0 ) ~ Cj Cf^F j(X k , s^, s^, Sjfc, Sjk) (14) 

k~l 

where k is the quadrature point index, n * is the number of quadrature points, Cfc is the quadrature constant associ- 
ated with point A;, X& = X(t*), s ik = s i(tk) and so on. 


We have imposed many types of constraints in the continuous time system. Due to space limitations, we’ll limit 
the discussion here to the most difficult of these constraints: nonlinear path constraints. Nonlinear point constraints, 
linear, and bound constraints are much easier to handle, and we leave them as an exercise for the reader. The 
nonlinear path constraints are constraints that must be satisfied along a section of the trajectory and not at a 
specific point. We can group the equality and inequality path constraints together for convenience as 


G(X 0 ) 


G £ (X 0 ) 

G/(X 0 ) 


/ Fifc(Xfc, Sifc, Sifc, Sifc, Sifc) \ 

[ F2fc(Xfc, Sjfc, Sifc, Sjfc, Sik) 1 


^^(X/e, S ik, S ik, Sik, $ik)) 


\ (Xfc , Sjfe, Sik, ^i/e)) / 


(15) 


where the function P ^ is the £ th path constraint function evaluated at the k th quadrature point. Hence, we now 
have (N ■ k) point constraints in the discrete time system, whereas we only had N path constraints in the continuous 
time system. 


NLP algorithms are almost always more efficient when they are supplied with analytic derivatives as opposed to 
using finite differencing to find approximations for the derivatives. Recall that the parameters we allow the NLP 
algorithm to vary are the initial cartesian states, X 0 . We need to determine the derivative of the cost and constraints 
with respect to X G . We can write the derivative of the cost function as: 


~ Cj Cfc^-Fj(X fe , Sjfc, s ik , Sik, Sik) 


k=i 


■ dX 0 


Using the chain rule we can expand Eq.(16) to read 


V ( dFj dXk dFj dSlk dFj dSlk dFj dSik dFj dilk ) 

dx 0 a Ck Ux fc ax 0 ds ik dX a ds ik dX a ds ik 0X o ds lk dx a ) 

k—1 ' 


(16) 


(17) 
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Similarly, for the constraints we can write 


<9G 

dX 0 


( • v \ 

( X k , S iky S ik , Sik ) Sik) 


dX 0 

dP 2k 

0X o 

dPek 

0X o 


{X k , S ikf S ik , Sik t &zk) 

( Xk , Sjfc, Sifcj Sik 5 


( 18 ) 


I 

V ^X G 


(Xfc, Sjfcj 


Sifcj Sjfc, Si/c)) 


/ 


Using the chain rule we can write the derivative of the I th constraint at the k th quadrature point as 

dP ik _ dPifc dX k dPadsik dPj k ds ik dP ek ds ik dP ek ds ik 

dXo ~ dx k dx 0 + ds ik ax 0 + ds ik dx Q + ds ik dx 0 + ds lk dx 0 


Inspecting Eqs. (17) and (19) we see that two types of terms appear and that these terms appear as pairs in a 
product. The first type of term is the derivative of the user supplied function with respect to the formation geometry. 
Examples of this type of constraint are 

dF J dFj &P g / 2n s 

8X k ’ ds ik ’ ds ik [ ) 

The second type of term is the derivative of the geometry at a specific quadrature point, with respect to the formation 
initial conditions. Examples of this type of term are 


dX k ds i k dsi k . . 

0X? 9X 0 ' dXo K } 

What makes this method “simple” is two- fold. First, the terms of type one contain derivatives of the mission 
performance metric function with respect to the formation geometry. For many missions, the performance metric 
can be cast as an algebraic function of the formation geometry. As we’ll see in the examples presented in later sections, 
these derivatives are trivial. The second reason this method is “simple” , is that the second type of terms contain 
derivatives of the geometry with respect to the initial conditions, and the form of these derivatives is independent 
of the mission specific performance metric, and therefore they only have to be derived and implemented once! It 
is very convenient that the derivatives are separated in this way, and the fact that they are is simply due to the 
application of the chain rule to a cost function that is explicitly dependent upon the formation geometry. Let’s look 
at the derivatives of type two, which are the derivatives of the formation geometry and its evolution, with respect to 
the formation initial conditions. 


Derivatives of Formation Relative Geometry with Respect to Orbit Initial 
Conditions 


The derivatives of the formation geometry with respect to the formation initial conditions contain information on 
the sensitivity of the state of the formation at some time, t k , with respect to the state of the formation at the 
initial epoch. We would expect portions of the STM to appear in these derivatives. Recall that the definitions for 
the geometry given in Eqs. (4-7) express the relative geometry, such as the relative position vector between two 
spacecraft, in terms of the cartesian states of the spacecraft. To calculate the derivative of side vector with 
respect to the initial position of the p th spacecraft, r op , we can write 


ds ik d , 

dr op dr op ^ jk * ^ 


A j(t k ,t 0 ) p^j 

—A t(t k ,to) p = £ 

0 otherwise 


(22) 


where A j is the upper left 3x3 component of the STM for the j th spacecraft and so on. 
subcomponents of the STM as 


*(Mo) 


A3;c3 (t, £ 0 ) ^3x3 (^5 ^o) A 

to) T>3x3 (Mo) J 


We define the four 3x3 
(23) 
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The derivative of with respect to v op can be written as 


B j(t k ,t 0 ) p = j 
-B e(tk,t 0 ) p = t 
0 otherwise 


ds ik 


d 


dv op dv op 


i r jk r (k) — 


The remaining derivatives of the formation relative geometry with respect to the initial condtions are: 

d s ik 


dr, 


Op 


dv OI) 


d 

dr 0 z 
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Cj{t k ,t 0 ) p = j 
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0 otherwise 
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Dj(t/t,t 0 ) p = j 
-D e (tk,t 0 ) p = i 
0 otherwise 

1/2 sL ds ik 
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Sik dr, 
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Sik dv, 
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d*ik 
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dsik , sf. dsik 


dv a 
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Sik dv, 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


op 


*op Sik \ *iks 

The derivatives contained in Eqs. (22-30) only contain terms that involve the spacecraft states and their associated 
STMs. Therefore, these equations do not change from one problem to the next. What will change is the actual 
numbers in the STMs, and the spacecraft states. This enables us to implement the algorithm in a general way in 
order to minimize the amount of recoding to apply the algorithm to a new mission or mission phase. Now let’s take 
a look at some aspects of the implementation we use. Then, we’ll look at two mission applications. 


Impementation 

The implementation of a numerical optimization approach must be carefully considered to ensure speed and ease 
of use. There are several aspects of the implementation that will have an influence on the speed of the approach 
we present in this paper. The first is how the orbit states and STMs are propagated. It is recommended that all 
propagation is performed in a compiled language as opposed to an interpreted language such as MATLAB. This is 
primarily due to the fact that we must propagate 42 differential equations per spacecraft. For our implementation, 
we coded most of the algorithm in MATLAB, and use MATLAB ’s fmincon SQP routine. However, all propagation 
is performed in compiled C code that can be called directly from MATLAB through a MEX interface. 

The second issue that influences speed is the number of quadrature points chosen, and the considerable amount 
of bookeeping which must be performed if the orbits are propagated and stored before evaluating the cost and 
constraints. For the examples presented in the next few sections, we used a simple trapezoid rule as opposed to a 
complex quadrature rule, in the summation of the discrete cost function. The number of points in the summation of 
the cost function was chosen to maximize the speed while yielding an acceptable approximation to the exact, integral 
form of the cost function. 

The second issue to consider is the organization of the code. To enable convenient application to new problems, it 
is helpful to isolate code that is specific to the user-defined function. By separating and organizing the code carefully, 
it is relatively simple to provide a new function containing cost and derivatives terms. For example, for the MMS 
example in the next section, there are only 16 lines of code that are specific to the MMS problem and they are all 
located in a single matlab function that is easy to modify to solve a new problem. 

Now that we have discussed a few of the important issues in implementation, let’s take a look at a specific mission 
example, the Magnetosphere MultiScale mission (MMS). 
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Application: MMS 

In the preceding portions of this paper we posed a general parameter optimization approach, and an approach to its 
implementation, to find optimal initial conditions of spacecraft formations. In this section, we’ll discuss an applica- 
tion of the method to the MMS mission. We’ll begin with a brief overview of the MMS mission. Then we’ll present 
the cost and constraint functions we chose for MMS and derive the necessary derivatives. We conclude the section 
with a discussion of results for a particular phase of the MMS mission. 

The MMS mission is one of several missions in NASA’s Solar Terrestrial Probes (STP) program. MMS will employ 
a four spacecraft formation to make fundamental advancements in our understanding of the Earth’s magnetosphere 
and its dynamic interaction with the solar wind. MMS will not be the first mission to use multiple spacecraft to 
study magnetospheric dynamics. The Cluster II mission was successfully launched by the European Space Agency 
(ESA) in the summer of 2000 and has already provided fascinating results on magnetospheric dynamics. There are 
several phases in the MMS mission. In this work we focus on Phase I, which consists of a highly elliptic reference 
orbit with a perigee radius of 1.2 Re and an apogee radius of 12 Re. The mission design objective is to provide a near 
regular tetrahedron formation between true anomalies of 160° and 200° . There are many ways to pose a set of cost 
and constraint functions to achieve this goal and in the next section we’ll discuss the method we choose. There have 
been many important contributions to the literature on tetrahedron formations and magnetospheric missions, and a 
detailed discussion of all the past efforts is beyond the scope of this work. References [1-7,9-14] contain information 
on the science of the MMS mission, reference orbit design, launch window analysis, tetrahedron design, and maneuver 
design among other topics. 

MMS Cost and Constraint Functions 

In Phase I of the MMS mission, one of the mission design objectives is to provide a near regular tetrahedron, with 
sides of 10 km, within a region defined by ±20° in true anomaly about orbit apogee. This requires casting a set 
of cost and constraint functions that take into account the size of the tetrahedron, and its likeness to a regular 
tetrahedron. Let’s begin by investigating some useful properties of tetrahedrons. If we assume one of the spacecraft 
is the reference, then there are three relative position vectors that describe the relative geometry of the spacecraft 
and we’ll define them as Si*, S2fc> and S3k to be consistent with our previous definitions. Knowing these quantities, 
we can calculate the volume of the tetrahedron using 

Va = ” |sifc ■ (S 2 k X S 3fe )| (31) 

where V a stands for the volume of the actual tetrahedron formed by the spacecraft, as opposed to the desired volume 
which we discuss below. Their are six unique sides in the formation, and the side lengths can be used to calculate 
the average side length, L*, using 


L* — - (sifc + S 2 k + «S3fc + $4k + $5k + S6k) (32) 

at time t*. Paschmann 10 presents several methods for evaluating how close a tetrahedron is to being a regular 
tetrahedron. One approach is to compare the volume of the actual tetrahedron V a , with the volume of a regular 
tetrahedron, V r , with side lengths equal to the average side length of the actual tetrhedron. The volume of a regular 
tetrahedron of side length L can be calculated using 

Vr = (33) 


Using V a and V r , Paschmann suggests an instantaneous volumetric performance metric, 
spacecraft as 





for a tetrahedron of 
(34) 


where the superscripts “fc” indicate values at the k th point. Q v has the useful property: 0 < Q v < 1. However, it 
does not take into account the actual size of the tetrahedron. We propose a simple polynomial function, Q s , that 
has the properties below, to take into account the size of a tetrahedron. The constants £i, £ 2 , ^3, and £4 are used to 
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change the shape of the funciton. A plot of Q s for a 10 km tetrhedron, with = 4, I 2 = 6, £3 = 18, and £4 = 25, is 
shown in Fig.l. The important property of Q s , is that it is zero for tetrahedrons with an average side length of less 
than 4 km or greater than 25 km. Between the values of 6 km and 18 km, Q s is equal to one. 


Ql(L*) 


' 0 

(L*-e 1 ) 2 (L* + c l -2e 2 ) 2 /(e2-h) i 

< 1 

(L* - £ 4 ) 2 (L* - 2£ 3 + U) 2 /{U ~ £s) 4 
K 0 


L* <£, 
li<L* < £ 2 
(2 < L* < £ 3 
£3 < L* < £4 
L" >£4 


(35) 


By making a composite function involving both Q v and Qs we can create a performance metric that evaluates both 



Figure 1: Plot of Q S (L*) 


the size and shape of a tetrahedron. Let’s use the following continuous-time form 

P(t) = 1 - Q v (t)Qs(t) (36) 

P(t) will be zero for a regular tetrahedron with a side length between 6 and 18 km. For unacceptable tetrahedrons, 
P(t) is one. We can formulate a cost function by integrating P(t) over the region of interest, or 

j= [ f P(t)dt= [ f (l - Q v (t)Q s (t))dt 
Jta Jto 

where t Q is the time when a selected spacecraft is at a true anomaly of 160°, and tj is the time when the reference 
spacecraft is at a true anomaly of 200°. However, we can’t evaluate this integral exactly, so we have to approximate 
it. We use a simple trapezoid rule, as opposed to a complicated quadrature rule, where 

J * c£(l - tfoi) A<* = C £(1 - 4 0,0*1^ (37) 


For MMS, the independent variables are the initial cartesian states of the four spacecraft in formation. In addition 
to minimizing the cost function above, we also must satisfy the nonlinear constraints that the semimajor axes of all 
spacecraft in formation are equal. We can write this as 


— 


1 op 


op 


(38) 


where a p is the semimajor axis of the p th spacecraft, ad is the desired semimajor axis of all spacecraft, r op is the 
magnitude of the initial position vector of the p th spacecraft, and v op is the magnitude of the initial velocity vector 
of the p th spacecraft. 
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In Eq. (37), we have a performance metric for MMS that is an explicit function of the relative geometry of the 
formation. In Eqs. (38), we have four nonlinear constraint functions in terms of the initial spacecraft states. Now all 
we require are the derivatives of the cost and constraint functions in order to apply the method to find an optimal 
solution. 


Derivatives of MMS Cost and Constraint Functions 

The last step to perform before solving for an optimal formation for MMS, is to determine the derivatives of the cost 
and constraints with respect to X G , s^, Sik, and The k th term in the summation in the cost function, J*., 
has the form 

Jk = (l - Q k MQ k s (sik )) (39) 

Notice that the cost function does not depend on the side rates, only on the side vectors and side lengths. Therefore 


^4 

dkij 


= 0 


and 


dJ fc 

dsn 


= 0 


The derivatives of the cost function with respect to s ^ are 


dJk _ ^ 2 Sifc • S2A; X S3fc 

dSi k |SU S 2 fc X S 3 fc| 


(8f*S z k fQ k 


dJk 

ds 2 k 
dJk 
dszk 


2 Sifc • S2fc X S 3 fc 

Q 

V2 L 3 fsifc S 2fe X S3* I 
2 Sifc • S2fc X s 3 fc 

— — -■ ■■ . ■ I 

y/2 L 3 |sifc • s 2 fc x Slfcl 


( s lfc s 3fc) Qs 
' ( S lfc S 2fc) Qs 


Now let’s look at the derivatives of the cost function with respect to s,*. We start with 


(40) 

(41) 

(42) 

(43) 

(44) 


dJ k 

ds ik 


“Si* 1 - 


= - O k ^ J - 

QvQs) Q^dsik 


dQt 

dsik 


Q 


k 

s 


We can show that 

dJk _ Qak Jsifc • s 2 fc x s 3 fc] 

ds ik 6 dL y/2 L *“ 

where 


dQl 

dL* 


0 

L* <h 

4(L* - €i)(L* + l x - 2t 2 ){L* - 4)/(~4 + h) 4 

CM 

V 

V 

0 

i 2 <l* < 4 

4(L* - 24 + U){L* - U){L* - 4)/(4 - 4) 4 

4 < L* < 4 

0 

A 

# 


Finally, the nonzero derivatives of the semimajor axis constraints are 


(45) 

(46) 


(47) 


da p 

da p dr op _ 0 / 

f a ' 

\ 2 r T 

\ l op 

dr op 

dr op r 0 p ^ 

\ r op > 

) Top 

da v 

dcip dv op 

2 Op 



0 V op dv op V 0 p fl 


(48) 

(49) 


We now have all of the information to use an NLP code to find and optimal solutions. In the next section, we 
discuss the results for an optimal MMS formation. 
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Results of Optimization for MMS 

Recall that the primary formation flying goal for the phase of MMS we consider is to provide a near regular tetrahe- 
dron, with side lengths near 10 km, in the orbit region between 160° and 200° true anomaly. We also must satisfy 
the constraint that the semimajor axes for all spacecraft are the same. Figure 2 shows several plots that illustrate an 
optimal solution. The top figure shows the evolution of Q s > Q v , and Q s ■ Q v for one Keplerian orbit. The lower plot 
shows the evolution of the six individual sides and the average side length. We see that the quality factor is near 
one during the region of interest, which is bounded by the vertical dashed lines. The optimal solution is significantly 
better that the initial guess which was calculated using an algorithm described in Ref [ 9 ]. For the initial guess, 
the time averaged quality factor Q v ■ Q s was about 0.55. The NLP algorithm found an optimal solution where the 
time- averaged quality factor was about 0.94. For this solution, the run time was about 35 seconds, and required 
about 40 function evaluations. From inspection of the bottom plot in Fig. 2 we see that the average side length is 
on the order of 10 km. The orbit states for this solution are found in the Appendix. 


Figure 2: Results of MMS Orbit Optimization 


Quality Factor vs. Elapsed Time in Orbit 





Now let’s take a look at a more complicated example, the Laser Interferometer Space Antenna (LISA), which has 
a large number of nonlinear path constraints. 
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Application:LISA 

LISA is a NASA/ESA mission to detect and study gravitational waves from massive black hole systems and galactic 
binaries. Gravitational waves are ripples in space-time caused by massive objects. They are a prediction of general 
relativity, but are yet to be directly detected. The detection and understanding of gravitational waves can provide 
breakthroughs and refinements in current relativistic theory. For a more detailed discussion of the LISA mission, see 
Refs. [8, 15-17]. 

The nominal LISA formation consists of three spacecraft in heliocentric orbits trailing the Earth by about 20°, 
with inclinations near 1° with respect to the ecliptic plane. The mission design goals for LISA are challenging. The 
primary goal is to provide a formation that maintains a nearly equilateral triangle with sides near 5xl0 6 km for 
the entire life of the mission, which is currently about 8 years. This has to be achieved entirely through careful 
orbit design, as continuous feedback control of the orbits is not permitted because it will interfere with the science 
measurements. We also must ensure that the sides of the triangle remain within one percent of 5xl0 6 km and that 
the side rates never exceed 15 m/s. There is a secondary, and competing goal, that we keep the formation as close to 
Earth as possible for power reasons. Let’s look at how we can cast these goals and constraints in a manner consistent 
with the algorithm presented in previous sections. 

LISA Cost and Constraint Functions 

There are essentially two goals for LISA that can be formulated as a cost function. The first is the desire to keep the 
side lengths as close to 5xl0 6 km as possible. The second is to minimize the distance between the LISA formation 
and the Earth. We choose to pose the cost function in terms of the Earth distance, and cast the side length goals as 
constraints. The following cost function permits us to minimize the distance of the LISA formation from Earth: 

n k 

J = c^2 (n k + r 2 k + r 3 k ) (50) 

fc=l 


where rik = (x^ + x\ h + and so on. To ensure that the sides of the formation never exceed a one percent 

variation from the desired side length, we can apply the following set of i ■ k path constraints. 

{sik ™ L ) 2 < gi{t k ) 2 (51) 


where L = 5xl0 6 and 



-20004 + 50, OOOfcm 
40, 000 km 


tk < 5 yrs 
t k > 5 yrs 


(52) 


where t k is the time at the k th point, in units of years. The function gi(t k ) above, and 52 ( 4 ) below, were suggested 
by Sweetser. 17 #i(4) and 52 ( 4 ) are chosen so that the range and range rate requirements are still satisfied with 
expected orbit insertion errors. The following set of i * k constraints ensures that the side rates never exceed 15 m/s: 


(sik) 2 < 52 (4) 2 


where 


92{t) = 


--tk + 15m/s 4 < 5 yrs 

5 

12 m/s 4 > 5 yrs 


Finally, we require that the ecliptic inclination, i ec , of each orbit be less than one degree. 


(53) 

(54) 


iec < 1 ° 


Let’s reformulate this constraint so that the derivatives are simpler. 

cos i ec = z T h ec > cos 1° 


where 

z = [0 0 1] T 


(55) 

(56) 

(57) 
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free — *"ec * Wc (58) 

and r ec and v ec are the initial position and velocity of a LISA spacecraft, in the Sun-centered mean ecliptic of J2000 
frame, and are given by 

r ec — *E/S + Rr 0 *;; (59) 

Vec = V E/5 + Rv o/c ; (60) 

where r E /s and 'Ve/s are the position and velocity of the Earth with respect to the Sun, at the initial epoch, and R 
is the rotation matrix from the Earth’s mean Equator of J2000 to the Mean ecliptic of J2000. 

The derivatives of the cost and constraint functions are presented, without derivation, in Table 1. The leftmost 
column in the table indicates the function, and the rows are the derivatives with respect to the term that appears 
in the column labels. The superscript “x” indicates the skew symmetric matrix. For derivatives with respect to the 
initial cartesian states, the “p” subscript is a dummy variable associated with the state of the p th spacecraft. 


Table 1: Derivatives of LISA Cost and Constraint Functions 


Function 


&Sik 

dr op 

dv op 

n k 

dc Y, (r lk + r 2 k + r 3k ) 

k= 1 

d(s ik - L ) 2 

d{s ik ) 2 

0 

2 (s ik - L) 
0 

0 

0 

2 Sik 

dr op ~ Ap{tk ' to) 

0 

0 

P± = B p (t k ,t 0 ) 

OY 0 p 

0 

0 

d{x T )h ec 

0 

0 

zT (~ V S/S R - Rv °p) 

z T (r* /s R + Rr *,) 


Let’s discuss some practical considerations before moving on to the results of the LISA optimization. The current 
launch date for LISA is Jan. 2015. The initial epoch we use is 01 Jan. 2015 12:00:00.000 UTC. The force model 
used in the LISA optimization includes point masses from all planets and the Earth’s moon. The coordinate system 
is Earth’s mean J2000 equator. However, the converged solution is presented in the Appendix in the sun-centered 
mean ecliptic of J2000. Now let’s look at details of a representative initial guess and optimal solution for LISA. 

Results of Optimization for LISA 

Figure 3 contains plots illustrating characteristics of the initial guess and converged solution for an 8.5 year LISA 
trajectory. The plots associated with the initial guess are in the left column, the plots associated with the converged 
solution are in the right column. The independent variable for each plot is the elapsed mission time in years. The 
first row of plots shows the side lengths between spacecraft, or s*. The second row of plots shows the side rates, or 
Si. Finally, the third row contains plots illustrating the evolution of the spacecraft-to-Earth distance. 

Let’s begin by taking a closer look at some properties of the initial guess. The first point we see is that the initial 
guess is not a feasible solution. We see that the maximum variation in range between any two LISA spacecraft is 
around 700,000 km. This is over an order of magnitude larger than the desired maximum range variation of 50,000 
km. Similarly, the maximum range rate for the initial guess is around 140 m/s, which is about an order of magnitude 
higher than the 15 m/s constraint. From inspection of the plots, we see that the initial guess is unstable. It turns 
out that the instability is due to the fact that the initial guess has been placed too close to the Earth. We see that 
the distance from Earth to the formation varies from a minimum of about 31x106 km, to about 53x106 km. 

The optimal solution is found to have significantly improved properties over the initial guess. Table 2 contains 
a comparison of relevant statistics for the two cases. We see that the maximum variation in range between any two 
spacecraft is 50,000 km for the converged solution. This indicates that the range constraint is active at the solution. 
The range rate between any two spacecraft is 13.7 m/s. The penalty we have paid to satisfy these constraints, is that 
the formation was moved farther away from Earth. The minimum and maximum distance from any LISA spacecraft 
to the Earth, is about 49x106 km and 64x106 km respectively. 

In summary, we have seen that the method we present here can solve a complex mission design problem, with 
multiple real-world constraints, in a deep-space flight regime. We began with a set of goals and requirements for 
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Table 2: Statistics for LISA Initial Guess and Converged Solution 


Trajectory 

max(\sik\) (m/s) 

max(sik) (10 e km) 

min(sik) (10* km) 

max(rik) (10 5 km) 

Initial Guess 

140 

5.7 

4.4 

53 

Converged Solution 

13.7 

5.04 

4.95 

64 


the LISA mission, and cast them mathematically as a set of cost and constraint functions. Next, we presented the 
derivatives of the cost and constraint functions with respect to the formation geometry. Comparing the initial guess 
to the converged solution, indicates that we can satisfy the current mission requirements for LISA, and that we don’t 
need a feasible guess to find a solution. Now let’s look at some general conclusions we can draw from this work. 


Figure 3: Initial Guess and Results for LISA Optimization 




Elapsed Time, yrs 




Elapsed Time, yrs 
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Conclusions 


We began this work with a long list of objectives that is motivated by the need for a practical, efficient method to 
solve complex, real-world guidance problems to support a diverse set of formation flying missions. The approach we 
developed and presented meets many, if not all of our initial goals. The approach is non-linear and not restricted to 
specific ight regimes or small inter-spacecraft separations. By carefully implementing the method, we can minimize 
the amount of work that is required to solve new problems. For example, the entire MATLAB implementation of the 
cost and constraint functions for the MMS example, neglecting code for propagation and the NLP solver, consists 
of only 354 lines of code. Only 36 lines are specific to the MMS problem. The reason so few lines of code are 
specific to the mission problem is due to the fact that we assumed the cost function is an explicit function of the 
spacecraft relative positions, velocities, ranges, and range rates. This simple assumption allows us to perform half of 
the derivation of the derivatives a-priori in a general way, and implement them in software only once. This reduces 
the work that must be performed to solve new analysis problems. Furthermore, the derivatives that can be calculated 
a-priori are the derivatives of the dynamics with respect to the initial conditions. The derivatives that are problem 
dependent, are often simply derivatives of algebraic functions with respect to the formation geometry. 

We illustrated the applicability of the approach to two different missions: MMS and LISA. Both of these missions 
have a complex set of cost and constraint functions that are explicitly dependent upon the formation geometry. We 
presented optimal orbit solutions for both missions. 
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Appendix: Optimal Orbit States for MMS and LISA 


Table 3: State Information for a Representative, Phase I, 10 km Tetrahedron Formation 


Property 

MMS1 

MMS2 

MMS3 

MMS4 

a (km) 

42095.7 

42095.7 

42095.7 

42095.7 

e 

0.81818 

0.81798411 

0.81799342 

0.81800317 

i (deg.) 

27.8 

27.800078 

27.801283 

27.798145 

w (deg.) 

15.000001 

15.008968 

15.002286 

14.982171 

fi (deg.) 

0 

359.99962 

359.98943 

0.014249285 

v (deg.) 

180 

179.99657 

180.00303 

180.00225 


Table 4: Heliocentric Mean Ecliptic J2000 , 01 Jan 2015 12:00:00.000 


OE 

LISA1 

LISA2 

LISA3 

LISA-Center 

a 

149457958 

149457278 

149457372 

149457536 

e 

0.018265676 

0.00949038651 

0.00887263317 

0.0086422348 

i 

0.890058699 

0.993559828 

0.993293398 

0.0663052608 

u 

90.6009132 

32.5150215 

148.24379 

272.489236 

n 

3.64713792 

126.757235 

240.28666 

182.897519 

V 

343.945562 

278.418085 

51.016237 

-16.9098743 
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